import numpy as np
import pandas as pd

np.random.seed(2)

#三角分解（或LR分解，或Doolittle分解）A=LU
def LU_decomposition(A):
    n = len(A[0])
    L = np.zeros([n, n])
    U = np.zeros([n, n])
    for i in range(n):
        L[i][i] = 1
        if i == 0:
            U[0][0] = A[0][0]
            for j in range(1, n):
                U[0][j] = A[0][j]
                L[j][0] = A[j][0] / U[0][0]
        else:
            for j in range(i, n):  # U
                temp = 0
                for k in range(0, i):
                    temp = temp + L[i][k] * U[k][j]
                U[i][j] = A[i][j] - temp
            for j in range(i + 1, n):  # L
                temp = 0
                for k in range(0, i):
                    temp = temp + L[j][k] * U[k][i]
                L[j][i] = (A[j][i] - temp) / U[i][i]
    return L, U


if __name__ == '__main__':
    # A=np.random.randint(1,10,size=[3,3])  #注意A的顺序主子式大于零
    A = [[4,2,1], [2,2,0], [1,0,3]]  # 举一个例子
    print("原矩阵A：\n",A)

    L, U = LU_decomposition(A)
    print("单位下三角阵L:\n", L)
    print("上三角阵U或R:\n", U)